DEMODULATION AND PHASE ESTIMATION OF TWO-DIMENSIONAL 
PATTERNS 



Field of the Invention 

The present invention relates to two-dimensional code patterns and more 
particularly to the automatic demodulation of ftinge patterns. The present invention also 
relates to the automatic estimation of pattern and texture orientation of an image. The 
present invention further relates to phase demodulation of a sequence of fringe patterns 
with arbitrary and unknovm phase shifts between them. Such patterns typically occur in 
optical interferograms, but can also occur in interferograms produced by more general 
waves such as acoustic, synthetic aperture radar (SAR), ultrasound, and x-ray More 
general image sequences, such as video, can in some cases also be modelled as fringe 
pattern sequences. 

Background 

Methods exist for demodulating fidnge patterns. Examples of such patterns are 
illustrated in Figs. 5 and 14A. However, conventional methods of demodulating closed 
curve fringe patterns, including circular and elliptical fringe patterns, inevitably produce 
ambiguities. Ambiguities are tlaen typically resolved by including additional constraints. 
These additional constraints may include restrictions on the smoothness of curves within 
the Mnge patterns, 

In addition to ambiguities other defects and errors may also be produced by the 
preliminary demodulation process. An example of this can be found in the classic Fourier 
(Hilbert) transform method (FTM). When the FTM is applied to a whole image, artefacts 
are produced, the artefacts being related to discontinuities introduced by the half-plane 
filter used in Fourier space. The artefacts result in errors in phase estimation, and 
specifically in regions where an angle of the fringe pattern is close to perpendicular with 
the half-plane discontinuity line, 

A method adopted to overcome this problem is to choose the half-plane filter m 
Fourier space such that the discontinuity line does not cut through the signal in Founer 
space. 

However, the frequency components derived firom circular or elliptical fringe 
patterns are also "circular" and no discontinuity line could be defined which does not cut 
through the frequency signal. 
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Other methods based on local (small kernel) demodulation may avoid the errors 
of the FTM, but typically fail to disambiguate correctly when there are strong 
perturbations ra the underlying fringe pattern. Features such as bifurcations and fringe 
endings represent typically strong perturbations , 
5 Images are often characterised by a number of featiires including pattern and 

texture. A number of image processing operations depend upon the estimation of the 
orientation of these image features. 

Image processing/analysis applications where the estimation of orientation is 
important include: 
.,10 • edge orientation detection; 

• orientation input to steerable filters; 

- f • pattern analysis; 

i3 • texture analysis, where the local orientation of texture is estimated for 

characterisation purposes; and 
5 • orientation selective filtering, processing, and enhancement of images. 

One such area where the estimation of orientation is particularly important is the 
a demodulation of fringe patterns or equivalently, AM-FM modulated patterns. One area 
where such patterns exist is in fingerprint analysis. Other areas include the analysis of 
naturally fonned patterns such as an iris diaphragm, ripples in sand, and fluid convection 
20 in thin layers. The primary reason why it is useful to know the fringe orientation when 
demodulating such a jfringe pattern is because it allows for 1 -dimensional demodulators to 
be aligned correctly. 

Another area where the estimation of pattern orientation is used is in image re- 
sampling and enhancement. 
25 Known methods for estimating pattern orientation include gradient-based 

methods, A serious deficiency with such gradient-based methods is that at the ridges and 
valleys of images, the gradient has zero magnitude. Gradient methods are therefore 
unable to provide information about pattern direction in such regions. 

Intcrferometry refers to experiments involving the interference of two light 
30 waves that have suffered refraction, diffraction or reflection in the interferometer. Such 
experiments typically involve the testing of lenses, mirrors and other optical components. 
The principle parameter of interest in intcrferometry is the phase difference between 
interfering beams. 



533242US.doc 



-3- 



Present methods of demodulating sequences of phase-related fiinge-pattems 
generally rely on a priori knowledge of the phase-shifts between each of the individual 
patterns. When these relative phase-shifts are available, it is possible to estimate the 
spatial phase by using a generalised phase-shifting algorithm (PSA). Li c^ses where the 
5 relative phase-shifts are not a priori known, the estimated spatial phase will be incorrect 
unless special error-reducing PSAs are used or methods to estimate the actual phase-shift 
between patterns are used. 

Error-reducing PSAs are usefiil where the deviation of the actual phase-shift 
from the expected phase-shift is small (typically less than 0.1 radian) and deterministic. 
10 When the phase-shift deviation is larger and/or noa-detenninistic (ie essentially random), 
then other methods must be used. 

Methods have been developed to estimate phase shifts between phase-related 
fringe patterns. Certain of the methods, based on statistical (least squares) estimation, 
requires at least 5 fi-inge patterns to work. Methods based on fitting the phase shifts to an 
15 ellipse also require at least 5 fringe patterns. Methods based on the statistical properties 
of spatial phase histograms have been proposed, but are heavily dependent on certain 
restrictive, and often unrealistic criteria. More recently methods based upon the 
correlation between patterns have been proposed but are very sensitive to the number of 
full fringes in the frames. These correlation methods assume that the fringes are 
20 essentially linear in the frame, which means only trivial fringe patterns can be processed. 
Methods based on the maximum and minimum variation of individual pixels in different 
- patterns typically require at least 15 frames to operate effectively. 

Accordingly, current methods of demodulating sequences of phase-related 
fringe-patterns are unsatisfactory in many situations, especially for a small number of 
25 phase-related fringe-pattems, and many algorithms do not work in the case of three or 
four frames- 
Disclosure of the Invention 
It is an object of the present invention to substantially overcome, or at least 
ameliorate, one or more disadvantages of existing arrangements. 
3C' According to a first aspect of the invention, there is provided a method of 

demodulating a real two-dimensional pattern, the method comprising the steps of: 

estimating a quadrature rwo-dimensional pattern from said real two-dimensional 
pattern using a two-dimensional spiral phase filter; and 
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creating a demodulated image by combining said real two-dimensional pattern 
and said estimated quadrature two-dimensional pattern. 

According to another aspect of the invention, there is provided an apparatus for 
demodulating a real two-dimensional pattern, the apparatus comprising; 
5 means for estimating a quadrature two-dimensional pattern from said real two- 

dimensional pattern using a two-dimensional spiral phase filter; and 

means for creating a demodulated image by combining said real two-dimensional 
pattern and said estimated quadrature two-dimensional pattern. 

Accordmg to another aspect of the invention there is provided an apparatus for 
10 calculating a quadrature conjugate of a real two-dimensional code pattern, the apparatus 
comprising: 

a first spatial li^t modulator for modulating coherent hght to produce said real 
two-dimensional pattern; 

a first lens for Fourier transformation of said real two-dimensional pattern to 
1 5 produce a first Fourier transformed image; 

a second spatial light modulator or spiral phase plate for phase modulating said 
first Fourier transformed image to produce a phase modulated image; 

a second lens for Fourier transformation of said phase modulated image to 
produce a second Fourier transformed image; 
20 ■ a third spatial Ught modulator or spiral phase plate for phase modulating said 

second Fourier transformed image to produce said conjugate of said real two-dimensional 
- pattern; and 

an image sensor for capturing said conjugate of said real two-dimensional 

pattern. 

25 According to yet another aspect of the invention there is provided a method of 

estimating an orientation angle of a pattern in an image, said method comprising the steps 
of. 

appiymg a complex energy operator to the image to provide an energy encoded 

image; 

^0 detennining a phase component of said energy encoded image; and 

calculating said orientation angle from said phase component of said energy 
encoded image. 

According to yet another aspect of the invention there is provided an apparatus 
for estimating an orientation angle of a pattern in an image, said apparatus comprises. 



533242US.doc 



means for applying a complex energy operator to the image to provide an energy 
encoded image; 

means for determining a phase component of said energy encoded image; and 
means for calculating said orientation angle from said phase component of said 
energy encoded image. 

According to yet another aspect of the invention there is provided a method of 
estimating a spatial phase of fiinge pattern images m a sequence of fringe patterns, said 
method comprising the steps of: 

a) removing offsets from each of said fringe pattern images to obtain pure 
AMFM patterns; 

b) determining contingent analytic images corresponding to each of said 
AMFM patterns; 

c) determining phase differences from dependent pairs of said contingent 
analytic images; 

d) estimating phase shifts between pairs of said contingent analytic images: 

and 

e) estimating a spatial phase of said fringe pattem images. 

According to yet another aspect of the invention there is provided an apparatus for 
estimating relative phase shifts between fringe pattem images in a sequence of phase- 
related fringe patterns, said apparatus comprising: 

means for removing offsets from each of said fringe pattem images to obtain pure 
AMFM patterns; 

means for detennining contingent analytic images corresponding to each of said 
AMFM patterns; 

means for determining phase differences from dependent pairs of said contingent 
analytic images; and 

means for estimating phase shifts between pairs of said contingent analytic 
images, wherein said phase shifts between pairs of said contingent analytic images are 
said relative phase shifts between fringe pattem images. 

Brief Description of the Drawings 

A number of embodiments of the present invention will now be described with 
reference to the drawings and appendix, in which: 

Fig. 1 is an illustration of the local structure of a fringe pattem; 
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Fig. 2 is an illustration of the Fourier transform of the local fringe pattern of 

Fig. 1; 

Fig. 3 IS a flow diagram of a method of phase demodulation; 

Fig. 4 is a flow diagram of an alternative method of phase demodulation; 

Fig. 5 is an example of a fringe pattern image; 

Figs. 6 and 7 are the real and imaginary parts of an intermediate result in 
applying the method of Fig. 4 to the image of Fig. 5; 

Fig. 8 is an alternative estimate of fringe orientation of the fringe pattern image 
of Fig. 5; 

Figs. 9 A and 9B are the magnimde and phase components of a demodulated 
fringe pattern when the estimate of fringe orientation of the fringe pattern image as shown 
in Fig. 8 is used; 

Fig. 10 is a schematic block diagram of an optical spiral phase demodulator; 

Figs. 11 A and 1 IB are the magnitude and phase components respectively of an 
intermediate result in applying an energy operator to the image of Fig. 5; 

Figs, lie and 11 D are the magnitude and phase components respectively of 
another intermediate result in applying the energy operator to the image of Fig. 5; 

Figs. HE and IIF are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 5; 

Fig. 12A is an example of a fringe pattern image containing noise; 

Figs. 12B and 12C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 12 A; 

Figs. 12D and 12E are the magnitude and phase components respectively of 
applying a modified energy operator to the image of Fig. 12A; 

Fig. 13 A is another example of a fringe pattern image containing noise; 

Figs. 13B and 13C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig, 13 A; 

Figs. 13D and 13E are the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig. 13 A; 

Fig. I4A IS an example of a fhnge pattern image with amplitude and frequency 
modulation; 

Figs. 14B and 14C are the magnitude and phase components respectively of 
applying the energy operator to the image of Fig. 14 A; 
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Figs. 14D and 14E axe the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig. 14A; 
Fig. 1 5A is an example of a texture image; 

Figs. 15B and 15C are the magnitude and phase components cespectively of 
applying the energy operator to the im^age of Fig. 1 5 A; 

Figs. 15D and 15E are the magnitude and phase components respectively of 
applying the modified energy operator to the image of Fig, 15 A; 

Fig. 16 is a flow diagram of a method of estimating the pattern and texture 
orientation of images; 

Fig. 17 shows the Fourier domain coordinates; 

Figs. ISA and 18B show the phase shift parameters d„ for the case where N=5 
on a polar plane; 

Fig. 19 shows a flowchart of an automatic calibration method; 
Fig. 20 shows yet another an example of a fringe pattern; 

Fig. 21 A and 2 IB show the magnitude and phase components respectively of an 
analytic image of the example fringe pattem shown in Fig. 20; 

Fig. 22A shows a phase difference between dependent analytic images; 

Fig. 22B shows a histogram of a spread in values of the phase difference shown 
in Fig 22 A; 

Fig. 23 A shows a modulus of the phase difference between dependent analytic 
images shown in Fig. 22A; 

Fig. 23B shows a histogram of a spread in values of the modulus of the phase 
difference shown in Fig. 23B; 

Fig. 24 shows a polarity function calculated from the phase difference shown m 
Fig. 22A and the modulus of the phase difference shown in Fig. 23A; 

Figs. 25A and 25B show estimations of an amplitude modulation and a phase of 
the fringe pattem shown in Fig. 20; 

Fig. 26 shows an estimate of an offset of the fringe pattem shown in Fig. 20; 

Fig. 27 shows a simple weight function calculated from the estimated errors in 
the analytic image shown in Figs. 21 A and 218; and 

Fig. 28 shows a schematic block diagram of a general purpose computer upon 
which arrangements described can be practiced. 
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Detailed Description including Best Mode 

Consider a two-dimensional fringe pattern with intensity of the fringe pattern as 

foUows- 

fix,y) = a(x,y) + b(x,y),cos[2iKiuo{x-xo}+vo{y-yo}yz] ~ (1) 
5 Position coordinates (x,y) tnay be continuous for an analog pattern or discrete 

for digital patt^s. A background level is denoted by a(x,y) while an amplitude 
modulation term is denoted by b{x,y), both assumed to be slowly varying ftmctions. 
Phase term z(xQ,yo) presents the local offset of the overall phase term in the square 
brackets [] of Equation (1). 
10 Referring to Fig. 1 where a 2-dimensional fringe pattern or AM-FM pattern, 

located aroimd a point {xo,yo) and formed by a number of fringes 10 is shown. A small 
region of interest is defined by Q. Within this small region n, the spatial frequency 
components uo and vg are slowly varying so that they are effectively constant around the 
point {xo,yoj, making the carrier a linear function of x and y. A normal 12 to the local 
1 5 fringes 10 is at an angle J3^ to the x-axis. 

A noise term is omitted from Equation (1), but may be present in real fringe 
patterns due to the occurrence of blurring, non-linearities, quantisation errors, smudging, 
scratches, cuts, dust, etc. 

The intensity pattern f(x,y) is (generally) a rapidly varying function of position 
20 coordinates (x, y) . 

The Fourier transform (FT) of the above fringe pattern region Q is: 
j\f{x,y)ex^-2m{ux+vy)}ixdy 

" (2) 
=4i{"-v) + (i?n(«-Wo,v-vJexrfz;ir]+^n(" + "o=v+Vo)exi^-z;irI)e-'''^'^^*'^ 

wherein i = V— 1 

From Equation (2) it is clear that the Founer transform of the fringe pattern / 
25 has 3 lobe centres, the first being at the origin and another two at spectral coordinates 
(wo,vo) and (-mo,-vo). The last two lobe centres are mirror images of each other with 
regards to the origin. The effect of windowing the fringe pattern /over a small region Q 
IS to blur the otherwise sharp lobes. The central lobe An, representing the zero frequency 
component, can be removed, and will therefore be disregarded in the following analysis. 
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Methods that may be used for removal of the zero frequency compotients include using a 
suitably matched high pass filter. The remaining spectral lobes are illustrated in Fig. 2. 

The local spatial frequency components «o and vo of the fringe pattern/ satisfies 
the following Equations: 



wherein ^ is an angle of the normal to the fringes with the x-axis and is the spacing 
between fringes. 

The conventional Fourier transform method (FTM) of fringe analysis is based on 
the assumption that the angle /3o is confined to a small range for the entire fringe pattem, 

10 typically much less than rc radians. However, when considering fringe patterns with 
circular, elliptical carrier patterns or other closed curves, the angle ySfc covers the entire 
angle range from 0 to 2n radians. The FTM would therefore be deficient in solving fringe 
patterns created with closed curve carriers. The reason for this failure is related to the 
inability of the FTM to isolate the two (mirrored) lobes of the FT. As noted above, the 

15 frequency components derived from closed curve fringe patterns are also "closed curves". 
Therefore, referring again to Equation (2), both the second and third terms will have 
"closed curve" spectral components, causing a total lobe overlap. 

The preferred implementation removes directional discontinuities in the Fourier 
processing of the fringe pattem /by smoothly encoding the fringe direction by utilizing a 

20 filter with odd (180° rotational) symmetry, but (unlike the FTM) a smooth transition 
between adjacent frequencies. Also, the filter used in the preferred implementation does 
not aher the strength of different frequency components, and is therefore a unit magnitude 
or phase only filter. 



5 




(3) 



25 



A (complex) Fourier representation P(u,v) of the filter preferably has the 
following anti-symmetry qualities: 



>P(w, v)-exp[i A(«, v)]=-P(-w,-v)=-exp[z(A(M, v)^)] 



(4) 



A preferred filter chosen has a phase A(w,v) which is set equal to the polar angle 
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A{«,v) = <f> (5) 



w=gcos(^) "I 

v=^sinW \ ' (6) 

The resulting filter is termed a "spiral phase filter" because the filter phase X 
increases linearly as the polar angle ^ increases. The magnitude of the filter is constant 
(so the filter is a pure phase filter) and the phase has the form of a spiral or heUcal surface 
with just one turn. An equivalent form of the filter is simply: 

The Fourier transform of the fringe pattern after the central lobe has been 
removed is as follows: 

Applying the spiral phase filter Pfu.v) to Equation (8) changes the lobe i vmmetry to give: 

= {B^{u-u^,v- vo) exp[i;i^]exp[z>] - 5^ (u + , v + Vo ) expH;};] exp[r<i])e-^"^"^"'^°> 

Equation (9) can be inverse transformed, as is shown in Eqxiation (10), by 
implicitly assuming that the polar angle <^u,v) « ± /3b remains constant over each lobe. 
This is equivalent to assuming that the lobe subtends a small polar angle ^, which is a 
reasonable assumption for typical firinge patterns/ 

gix,y) = T r G(w,v)exp[2;n(wx + v>0]i^«'iv 

= z exp[i>9o ] sin[27r(M, {x- x^) ^v^iy - y^))^ x\ 
The desired sine component of the local fringe pattern gix,y') can then be 
extracted by multiplying si^,y) with an orientational phase component, exp[-//9b] as 
follows 

g(x. y/exp/ W^, ; = z sinp;r (uo{x- Xa}-rv,{y- y,})-^ xl ( H ) 
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Conceptually the process is repeated for all fringe regions with their 
corresponding fringe orientation emgles f}. In practice the procedure can be applied to all 
regions simultaneously. It is noted that the original (conceptual) windowing is not 
actually required. 

The aforementioned process depends upon the functions h(x.y), fi(x,y), z(^-y)> 
Uo(x.y), Vo(x,y), and <yo(x,y) having suitably slow variation. However, it is found in 
practice that the preferred implementation can be successfully implemented on any fringe 
pattern /with a suitably smooth variation in orientation fi(x,y) and spacing Ko. For fringe 
patterns / with circular fringes and constant or quadratic (chirp) fringe spacing Xo, the 
preferred implementation only breaks down near the very center of the fringe pattern / 
due to the assumption of small polar angle subtense failing near the centre. 

Equation (11) represents an estimation of the imaginary part of the analytic 
function (or image) related to the original real fringe pattern/ This component is also 
called the quadrature component or analytic conjugate. By forming a new (complex) 
output image comprising the original real image / and the estimated imaginary image, the 
phase and magnitude of the resulting complex fringe representation can be obtained as 
follows: 

<f- a)+{g{x,y) ex.p{-iPo]) « ^-exp{r[2;zi(«o{j:-xo}+vo{y->^o})+;^]} (12) 

Fig, 3 is a flow diagram of a phase demodulation process 1 00. The fringe pattern 
image f(x,y), captured in step 110, typically includes only a real part. The fringe pattern 
wnzgcf(x,y) is transformed in step 120 by applying a 2-dimensional FFT. Because the 
fringe pattern image f(x,y) is real, the real and imaginary parts of the transform F(u,v) are 
symmetric and antisymmetric respectively. 

The DC component (lobe) of tlie transform F(u,v) is removed in step 130. In 
step 40, the spiral phase fiherPfw.v) is applied to the output from step 130. In step 150, 
the inverse Fourier transform is calculated. This may be performed by again applying a 
2-dimensional FFT. 

In step 160, an oiientational phase filter is appUed to the resulting image of step 
150. In step 170 the real part is discarded, leaving only the imaginary part, which is then 
added to the input fringe pattern image f(x,y) in step 1 80 to produce an output image in 
step 190. 

Fig. 4 is a flow diagram of an alternative implementation of a phase 
demodulation process 200, and will be explained with reference to an example fringe 
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pattem. An example of such a two-dimensional fringe pattern is illustrated in Fig. 5, 
wherein only the real part of the pattern is known. 

The fringe pattern image f(x,y) is captured in step 210. In the phase 
demodulation process 200, steps 120 to 150 of the phase demodulation process 100 
shown in Fig. 3 are replaced by a single convolution step 220. The fringe pattern image 
f(x.y) is preferably convolved with spiral phase function p(x,y), wherein p(x,y) is defined 
as follows: 

Pi^.y) = _Q/'(tt,v)exp[2TO(MA:+ vy)]dudv 

= 1^ / 3 \ exp[2m(ux + vy)]ductv 



x = r cos 6'! 

>- = 7-sin^J ^^^^ 

The result after such a convolution on the fringe pattern image f(x.y) of Fig. 5 is 
illustrated in Figs. 6 and 7, showing the real and imaginary parts otf(x,y)*p(x,y). 

Referring again to Fig, 4, in step 230 the orientational phase filter, exp[-i3] is 
applied. However, an estimation of the fringe orientation p is needed. 

One method of estimating the fringe orientation p is to compute the 2-D energy 
operator which gives uniform estimates over the fringe pattern un^gc f(x,y) . In contrast, 
typical gradient based methods fail to give reliable orientation estimates at the peaks and 
valleys of the fidnges. 

An alternative method of estimating the fringe orientation |5 is to use an unifonn 
orientation estimator, but to utilise local enviroimiental constraints (such as topology of 
smoothly varying fringes) to resolve or unravel the orientational ambiguity. Fig. 8 
illustrates a smooth estimate of the fringe orientation exp(i^) derived from exp(2i^). The 
orientation phase estimate exp{ij3) has only one discontinuity line at positions where the 
phase changes from 360» to 0°. A number of different methods are suitable for resolving 
the orientational ambiguity. The magnitude and phase components of the result of step 
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230 using the estimation of the fringe orientation p of Fig. 8 is illustrated in Figs 9 A and 

9B. 

It is noted that an estimate of the orientation can be obtained from the output of 
step 150 in Fig. 3 or alternatively from the output of step 220 in Fig._.4. However, 
5 additional unravelling will be required to obtain a useable orientation estimate. 

A preferred implementation of estimating the fringe orientation p 4s described in 
detail below. 

In step 240 the real part is discarded, leaving only the imaginary part, which is 
then added to the input fringe pattern image f(x,y), in step 250 to produce an output image 
]0 in step 260. 

The method of Fig, 3 can be practiced using an optical spiral phase demodulator 
300, such as that shown in Fig. 10. The optical spiral phase demodulator 300 is 
illuminated by coherent light 301 . A first spatial light modulator (SLM) 302 modulates 
the coherent light 301 to produce the input image, which is typically a fringe pattern 
^ 5 f(x.y). The fringe pattern f(x,y) is Fourier transformed by lens 303. A second SLM 304 or 
a spiral phase plate performs phase modulation in the image projected on it. The phase 
modulated image is then Fourier transformed by lens 305 and is projected onto a fmal 
SLM 306 or spiral phase plate. The final spiral phase plate 306 retards the coherent 
optical beam in a manner such that the plate introduces a spiral phase multiplication. An 
20 image sensor 307 is in close proximity to the final SLM 306 to detect an output image. 
The output image will be in quadrature (or analytic conjugate) to the input image. In 
other words, if a cosine pattern is input, a sine pattern is output. 

In the special case where the input image is a circularly symmetric fringe pattern, 
then the SLM 306 has the simple form of a spiral phase with opposite sign to SLM 304. 



Returning to the description of the estimation of the fringe orientation p, ii is 
noted that in recent years, an efficient method of signal demodulation using a concept 
known as an "energy operator" has become popular. The operator can be used for 
demodulation of AM and FM signals, including combined AM-FM signals, using a result 
from the analysis of simple harmonic motion (SHM). The energy operator ^ is defined 



25 



as: 




(16) 
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A particularly useful property of the energy operator W is that it gives uniform 
estimates of modulation parameters anywhere on the signal / Essentially this is due to 
the fact that, at local peaks (and troughs) of the signal /, the energy is concentrated in the 
potential energy component, represented by the second term in Equation (i-6), whereas at 
the zero crossings, the energy is in the kinetic component represented by the first term m 
Equation (16). In other words, the rapidly varying sine and cosine components associated 
with AM-FM signals are exactly cancelled out by the energy operator ^, leaving an 
underlying constant representing total energy. 

The operation of the energy operator ^ has also been extended to N-dimensions. 
This is described in the publication "A multi-dimensional energy operator for image 
processing", SPIE Conference on Visual Communications and Image Processing, Boston, 
MA, [1992], 177-186, by P. Maragos, A-C, Bovik, and J. F. Quaxtieri. 

Utihsing the energy operator T on images, (i.e. data sets formed in 2- 
dimensions,) separate x and y components may be calculated and added together to obtain 
an overall 2-dimensional energy operator. 

However, the energy operator defined in this simple way has certain 
deficiencies. The main deficiency being that the orientation associated with the energy 
components is strictly limited to the positive quadrant of the x-y plane, which means that 
visually very different orientations are indistinguishable using two positive energy 
components. 

Referring again to Fig. 1 and Equation (1), the intensity of the fringe pattern can 
be represented as follows: 



wherein any slowly varying (non-multiplicative) backgrotmd offsets a(x,y) are removed 
with an appropriate high-pass filter. The local spatial frequency components uo and vo 
satisfy Equations (3). 

The X and y component energy operators and SKy respectively may be 
performed on the fringe pattern defined m Equation (17) as follows: 



g(x,y) = b{x,y).cos[2n(uo{x-Xo}+\o{y-yQ})+X] 



(17) 




(18) 
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^M'))-[^] -^0 (19) 

Substituting Equations (18) and (19) into Equations (3), the local frequency and 
orientation of the local firings pattern may alternatively be expressed as follows: 



(20) 



tm{p,) = ±^^^ (21) 

Equation (21) indicates a loss in sign information, which corresponds to an 
inability to distinguish between fringes at four different orientations 
^0'- /^(» ^ ^0'^- ^0- This is very undesirable because the human eye can clearly 
separate orientation JS^ from 7z-p^, and -^^ from Tc + fi^. In general, however, it is 
not possible to distinguish (on a local level at least) between any orientation and the 
rotated orientation. 

An alternative approach to extending the energy operator *F to two dimensions is 
by defining a two dimensional complex operator which encodes the energy in the 
magnitude and the orientation in the phase: 

{5} + [g] = {in>J {^l + ivl ) (22) 

As mentioned above, the simple x and y components alone fail to represent the 
orientation correctly. A preferable approach would be to have a 2D energy operator that 
has the following output: 

{s{x)\ = {inbf {u, + ^v, J ^ {2Ka,bJ exp(2^7?„ ) (23) 

The complex energy operator enables the isolation of the fringe spacmg Xo, as 
seen in Fig. I, and orientation angle po respectively. It is noted that the term 2 in the 
exponential factor of Equation (23) again causes the orientation angle (3o to be defined 
modulo -K. Such an energy operator is developed by defining a 2- dimensional complex 
differential operator D as follows: 



8y 



(24) 



- 16- 

The complex energy operator Tc can then be defined as follows; 

{D{g]f -gD^g]^{2KCT,bf^^p[2iP,)=^^^{g] (25) 
With the Fourier transform defined as follows: 

G {u,v)= \ ^g{x,y)Ql^'p[-2m{ux + vy)^dy (26) 

the operator D can be implemented in the Foxixier space , providing the following 
relationship: 

(2m X« + iv)G {u, v)= J jZ>{^{x, y)]exv[- 2m{ux + yy)}ixdy (27) 

wherein 

27r{u + iv) = 27Eq. exp(z>) (28) 
and the polar coordinate transforms are: 

M = ^COS(9))1 

v^qsm.{<p) i (29) 

The corresponding Fourier operation of Equation (27) provides for the 
multiplication of the Fourier transform of the fiinge pattern G(u,v) with a spiral phase 
filter, g&xip{i(p) The spiral phase filter has the characteristic of linearly (conical) 
increasing in magnitude with fi-equency. The Fourier transformation pairs may be written 
as: 

G{u,v)<^ g{x,y) (30) 
{27d)q exp(i>)G(w, v) <^ D{g{x, y)] (3 1 ) 

{irrif exp(2z>)G(w, v) ^ {gix, y)) (32) 

The derivatives in Equation (25) can therefore be implemented as 2-dimeusional 
Fourier filters with spiral phase and conical magnitude. 

By applying the differential operators defined in Equations (30) to (32) as filters 
through simple multiphcation of the image of Fig. 5, the results as shown in Figs. 1 1 A to 
UP IS obtamed. Figs. IIA and IIB show the magnitude and phase components 
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respectively of the (D{g}f term of Equation (25), calculated by using the relationship in 
Equation (31). Similarly, Figs, llC and UD show the magnitude and phase components 
respectively of the gD^ {g} term of Equation (25), calculated by using the relationship in 
Equation (32). The magnitude and phase components respectively of the complex energy 
operator is shovm in Figs. HE and UF, Referring to Equation (23), the phase 
component of the complex energy operator has double the orientation angle fio of the 
fringe pattern / 

All the above calculations and derivations may alternatively be implemented as 
convolution operators without Fourier trans fomiation. 

However, when the energy operator is applied to a fringe pattern / containing 
noise, such as that shown in Fig. 12 A, the resulting magnitude and phase components of 
the complex energy operator shown in Figs. 12B and 12C, are less satisfactory. The 
primary reason for the degradation is the linear (conical) spectral enhancement of the 
differential operator D. The differential operator D effectively enhances high frequency 
noise, as it essentially operates as a high-pass filter. In particular, in the areas of low 
spatial frequency ctq, such as the area 20 on Fig. 12A, the contribution by the noise signal 
due to the differentiation dominates the determination of the orientation of the firinge 
pattern. Referring to Fig. 12C, the phase component, from which the orientation can be 
obtained by using the relationship jfrom Equation (23), provides near random noise in the 
corresponding area 22. 

One possible correction is to apply a low pass filter function to the input image 
or the result. However, the energy operator ^ is not a linear (shift invariant) operator It 
is dependent on products of the input function g and its derivatives 8g/dx and d^g/dx^, as 
can be seen from Equation (16). Low pass filtering results in the image becoming 
blurred. 

To overcome the deficiencies of the complex energy operator while 
maintaining the ease of calculation of the orientation angle ^o, a pure phase operator I>^, 
thus having uniform (spectral) magnitude, is created by setting q in Equations (31) and 
(35) to unity, to define the following Fourier transform pair: 

exp(r«3)G « Dj^ {g(x, ;;)) (33) 
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The pure phase operator is a spiral phase (or vortex) operator which results 
in the removal of noise without significant blurring of the signal. By appljdng the 
modified energy operator Tm defined as: 



to the image in Fig. 12A, the magnitude and phase components respectively of the 
modified energy operator Tm shown in Figs. 12D and 12E, are obtained. The comparison 
of the phase components of the complex energy operator shown in Fig. 12C and 
modified energy operator shown in Fig. 12E, and in particular to the corresponding 
areas 22 and 24 respectively, indicate that the modified energy operator provides a 
better result. 

It is noted that, similar to the energy operator the modified energy operator 
SKm is also non-linear. A consequence of this is that it is generally not possible to obtain 
the performance of the modified energy operator by merely applying a low pass filter 
to the input function, nor by filtering the output function, of the energy operator ^. 

The modified energy operator has a property of scale invariance, which it 
inherits ft-om the scale invariant Dm operator. Scale invariance is another aspect of the 
operator's spectral neutrality and means that the operator essentially gives outputs 
independent of feature size. 

To further illustrate the advantages of the modified energy operator ^m. the 
complex energy operator and the modified energy operator Tm are respectively 
applied to an image with uniform noise in the range ±64 greylevels shown in Fig. 13 A. 
The magnitude and phase components respectively of the complex energy operator are 
shown in Figs. 13B and 13C, while the magnitude and phase components respectively of 
the modified energy operator are shown in Figs. 13D and 13E. Because of the high 
level of noise in the input image, the phase component of the complex energy operator 
is completely weighed down by the contribution of the noise. The modified energy 
operator gives a much clearer phase without any obvious blurring that would be 
expected of noise reduction by conventional low pass filtering of the complex energy 
operator 

Another comparative example wherein the performance of the modified energy 
operator can be seen is shown in Figs. 14A to 14E, Fig. 14A shows an amplitude and 
frequency modulated image. The magnitude £ind phase components respectively of the 




(34) 
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complex energy operator are shown in Figs, 14B and 14C, while the magnitude and 
phase components respectively of the modified energy operator Tm are shown in Figs. 
14D and 14E. The complex energy operator again has difficulty in regions of low 
spatial frequency ?io, such as area 30, causing the corresponding area 32 in Fig. I4C to 
have high levels of uncertainty. On the other hand, the corresponding area 34 when the 
modified energy operator Tm is applied, provides a much clearer result. 

It is further noted that since the underlying fringe pattern of the images in Figs, 
12A, 13A and 14A are the same with different processes applied thereto, the processes 
not affecting the orientation of the fringes of the underlying pattern, the results of the 
modified energy operator shown in Figs. 12E, 13E and 14E are substantially similar. 

The modified energy operator may be implemented entirely in the Fourier 
space domain by using complex convolution kernels that correspond to the pure phase 
spirals. With the kernel k and its Fourier transform K, the Fourier transform pair can be 
defined as: 

Kiu,.) = exp(,>) . ^ <=> %tf) . i^P(i^ = y) (35) 
wherein the spatial polar coordinates satisfy the following: 

x = rcos(^) |- (36) 
y = rsin(^) 

The convolution kernel k reduces by the inverse square of the radius, and also 
has a spiral phase structure, but rotated by 90** with respect of the x-y coordinates. 
According to the above, since it is now known those operations required to be performed 
to extract local orientation information from a fringe pattern, it is not necessary to 
perform such calculations in the Fourier domain. Rather, the desired functions may be 
recalculated in the real domain thereby facilitating ease of processing. In many cases the 
desired functions may be approximated by suitably truncated kernels. The inverse square 
drop-off allowing relatively small truncation errors even for small kernels 

Referring again to Fig, 1, by considering a sufficiently small region Q of images 
composed of line or edge features, the local structure fits the model of Equation (1). This 
is true for many image features, the most obvious exceptions being texttues that contain 
niultiple AM-FM modulated signals. Sometimes an image may approximate the model of 
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Equation (17) except for a low or zero frequency background offset in intensity. In such 
cases this background level may be removed by high pass filtering (or related methods) 
leaving an image tliat does satisfy Equation (17) which can then be processed using the 
complex energy operator formahsm. 

Where Equation (17) is representative, the modified energy operator as 
defined in Equation (34), may be used for the calculation of the orientation -angle fio of all 
such images. 

An example of such an image and the operation of the modified energy operator 
are shown in Figs. 15A to 15E. Fig. 15A shows a texture image. The texture image 
is characterised by distinct lines 40 and areas of low variation 41. The magnitude and 
phase components respectively of the complex energy operator are shown in Figs. 15B 
and 15C, while the magnitude and phase components respectively of the modified energy 
operator are shown in Figs. 15D and 15E. The phase component of the modified 
energy operator SPm provides a result that is more intuitively representative of the 
orientation of the texture image of Fig. 15 A. 

A flow diagram of a method 400 of estimating the pattern and texture orientation 
of images is shown in Fig. 16. An input image is inserted in step 405. The input image is 
typically digital format with a predetermined number of columns and rows. Alternatively 
the image may be scanned to convert it into the digital format. 

Although the preceding analysis has been in terms of continuous fttnctions and 
continuous operations, the analysis is also directly appUcable to discretely sampled 
functions (images) and discrete operations. 

A Fourier transform of the image is calculated in step 410, typically using the 
Fast Fourier Transform (FFT). The Fourier transformed image is multiplied with a tenn 
exp('i<p) in step 420. 

An output of step 420 is again multiplied with the term exp(i(p) m step 430 
before it is inverse Fourier transformed, typically using an Inverse Fast Fourier Transform 
(IFFT) in step 440. The result is then multiplied with the input signal in step 470. 

The output fi-om step 420 is also inverse Fourier transformed in step 450 before it 
is squared in step 460, The output from step 470 is subtracted jSrom the output of step 460 
in step 480. From this complex signal, the angle of the signal is calculated in step 490 
before it is halved in step 494 to provide an orientation angle as an output in step 496. 
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As noted above, the orientation estimate is modulo n but can be imwrapped using 
continuity constraints to obtain an estimate that is modulo Ztc. 

Interferometry is the use of interference phenomena for measurement purposes. 
This allows for the measurement of very small angles or small displacements of objects 
relative to one another. An interferometer is a device to make such measurements, which 
typically uses a beam of light, coming from a single source, such as a star, a laser or a 
lamp, and two or more flat mirrors to split off different light beams. These beams are 
then combined so as to 'interfere' with each other. The beams interfere with each other 
by constructively adding together and cancelling each other out, thereby forming 
alternating bands of light and dark, called fringe patterns. 

What makes the interferometer such a precise measuring instrument is that these 
fhnges are only one light-wavelength apart, which for visible light corresponds to about 
590 nm. With the fringe patterns having a cyclic nature with phases that shifted, the 
important measurement parameter is the relative phase-shifts between these phase-related 
fringe patterns. 

Consider an intensity Mx,y) of the n'^ fringe pattern in a sequence of phase- 
related fringe patterns with the following mathematical form: 

fn {x, y) = y) + b{x, y)cos[zix>y)+ <5„ ] (37) 
wherein a{x,y) is again denotes a background or offset, b(x,y) the amplitude 
modulation, and Ji^x.y) the phase. Each fringe pattern has a phase shift parameter S„ . An 
example of such a fringe pattern is shown in Fig. 20. 

An improved method of estimating relative phase-shifts between phase-related 
fringe patterns f„ is proposed. The method is based on an isotropic quadrature transform. 
The isofropic quadrature fransform is used to estimate the spatial phase [jdx,y)+S„] in 
one fringe pattem/,, followed by comparing the spatial phase [zix,y)+3„] between fringe 
patterns /„ to estimate the relative phase-shift(s). With 3 or more frames, the method is 
independent of background variations a(^x,y) . 

The exact value of phase %(x,y) can be calculated from three or more fringe 
patterns;^ if the values of the phase shift parameter 5„ is known for each fringe pattern/,. 
For more than three fringe patterns f„ix,y) , the phase ;^x,y) is over-determined and is 
usually estimated in a maximum likelihood process. Other estimation processes may be 



553242US.doc 



-22- 



more appropriate in cases where additional systematic eixor or distortion effects apply to 
Equation (37). 

Consider a least squares estimate for the phase z(x,y) from N fringe patterns: 





(■Lf.) 











(38) 



where c„ =cos<5„, and =sin<5'„, and the summations are over N fringe patterns ie. 
^ - Least squares Equation (38) can be inverted, providing c„ and are such 

that a singular matrix does not occur, as follows; 
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(39) 



Equation (39) allows for the offset a(x,y), the amplitude modulation b{x,y), 
and the phase xi^.y) to be calculated. However, before Equation (39) can be used, all the 
phase shift parameters S,^ must be known. 

The phase shift parameter 5^ of each fringe pattern is estimated by using a 2-D 
quadrature technique based on a spiral phase Fourier operator. Estimation of phase xk^.y^ 
in the above manner has many advantages related to the extraordinary demodulation 
properties of the spiral phase Fourier operator. 

A first step in determining the phase shift parameters <5„ is to remove the 
background (offset) a{x,y') fix)m each of the fringe patterns fnix^y) . In the preferred 
implementation, this is done by forming inter-frame differences between pairs of fringe 
patterns fr,{x,y) as follows: 

= /, - /™ = = b{Q.o^ + S„\- co^\x + 5„ ]} (40) 
The inter-frame difference function g,,^ is a pure AMFM (amplitude modulation 
and frequency modulation) fimction without offset and therefore can be demodulated if 
the inter-frame difference function g^^ , which is also a fringe pattern, is locally simple. 
For the fringe pattern g^^ to be locally simple, it has to have a unique orientation 
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fi{x,y) and spacing ^{x,y) at each location {x,y), except, perhaps, at a finite number 
of locations where the orientation P{x,y) and/or spacing A,{x,y) may be uncertain. 
Advantageously, local simplicity is a property which distinguishes fringe patterns from 
other, more general, patterns. Fig. 1 shows the local structure of a typical fringe pattern 
for the small region Q. 

A spiral phase Fourier operator V{}, hereinafter referred to as a vortex operator, 
is defined as follows: 

V {/ {x, y)] = F-^ {exp {/ {x, y)}} 

F {/ y)} ^F{u,v)=jjf{x,y) exp [-2jrz {ux + vy)] dxdy 

F~'{F{u,v)] = f{x,y)=^ J I F(u,v)e1^p['^27ri(ux + ^y)']c^udv 

u = g cos <p} 

V = qsin<^ j 

The Founer operator F{ } transforms a general function /('x,>;; from the spatial 
domain {x,y) to the spatial frequency domain (w,v). The Inverse Fourier operator 
F"' { } transforms a spatial frequency domain signal back to the spatial domain {x,y). 
In practice, for discretely sampled images, the Fourier and the Inverse Fourier transforms 
are efficiently performed by a Fast Fourier transform algorithm and an Inverse Fast 
Fourier transform algorithm. 

The spiral phase exp[i^] is wholly a function of the angular component of the 
spatial frequency polar coordinates shown in Fig. 17. The vortex operator V{} is 

applied to the inter-frame difference functions g^^ to obtain: 

V{^„,„ } = i exp[zy0>[sin(;ir + Sj- sinijc + S„ )] (42) 

Prom Equation (42) it can be seen that the vortex operator V{} is an isotropic 
quadrature operator in that it converts cosine fringes at any orientation J5{x,y) into sine 
fringes at the same orientation J3{x,y). The approximation in Equation (42) is valid 
except in regions where the amplitude modulation b(x,y) , spacing X{x,y) or orientation 
P{x,y) vary too rapidly. Axeas where the radius of curvature of the fringe pattem is 
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smaller than the ftinge spacing X[x,y) also give a lower accuracy approximation. 
However, for most typical fidnge patterns g„„ , the approximation is good and any errors 
are isotropically distributed. In contrast, conventional methods, such as the half-plane 
Hilbert transform, cause highly directional errors, and are not suitable for'liigh precision 
5 demodulation of general ftinge patterns. 

Next, the local orientation /5 (x, y) is compensating for. To do so, the orientation 
^ [x, y) must be estimated. There are a number of ways to do this. One method relies on 
the gradient of the fringe pattern g„,„ . The ratio of the x and y components gives the 
tangent angle I3{x,y) from the local Taylor's series expansion of the phase xi^^y) about 
10 the point {xf^,yQ) : 

The gradient components are: 

^ . 6(- sinO, <5 J + sinU 5 J)^ 

^ (44) 
^ b{- sinC + S„ sinC + ))^ 

oy dy 

The ratio of the gradient components gives the tangent angle: 

^/^.^...,,. 

Again the approximation of Equation (45) is valid in regions with smoothly 
varying parameters as defined for the vortex operator V{}. Other methods of orientation 
estimation may be used. For example, methods based on the multi-dimension at energy 
operator or the vortex operator itself may be used. Generic orientation estimation is 
20 assumed at tliis point, and the orientation /3 for any fringe pattern /„ or inter-frame 
difference function may be calculated. Alternatively the orientation may be 
estimated from combinations of the previous methods in a statistically advantage manner. 
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It will be assumed hereinafter that there is just one orientation estimate, which shall be 
called orientation estimate p^. The orientation estimate y?^ typically contains an 
ambiguity in the sign: 

exp[z^] = cos yff + i sin = ±^j^ (46) 

Different orientation estimates typically have the sign ambiguities in 
different places. It is worth noting that the overall sign ambiguity arises from the general 
difficulty in distinguishing a fringe at orientation from a fringe at orientation . 

Equation (40) is rewritten as: 

g„„ = 2b &in[z + (S„ + S„ )/2]sin[(<5„, - )/2] (47) 

and Equation (42) as: 

exp[- iJ3, ]V{^„„, } = lib exp[/05 ~ y?, )]co^\;c + + S„ )/2]sin[(j„ - S„ )/2] (48) 

The complex exponential in Equation (48) has possible values {-1;1}, dependmg 
on the orientation estimate , and shall be called a polarity factor h{x,y)] 

h{x,y) = cxpli{/3-j3^)] (49) 

A contingent analytic image is defined as: 

- K^..„ - exp[- ip^ ]V{£_ }) (50) 

By substimting Equations (47), (48) and (49) into Equation (50), the contingent 
analytic image is rewritten as: 

i_ = 2b sm[(S„ -Sj/ 2lh cos[z + (S„ + <5„ )/ 2] + i sm[z + (S^ + S„ )/ ij) (51) 

This definition of an analytic image is contingent upon the orientation function 
used in the definition shown in Equation (50). It is noted that the contingent analytic 
image f ,„„ can be achieved by placing the inter-frame difference function g„^, in the 
imaginary part of a complex image, and placing the imaginary output of the vortex 
operator V{} in the real part. 

A phase a„„, of the contingent analytic image ^„,„ is calculated as follows. 

«^ = Arg{ig„„ - i exp[- ip^ JV{^„„ )]^h]j:+{Sm + Sn)/2\ (52) 
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The calculation of phase a„„, is possible as long as sinj 



- )/2] ^ 0 , which 



only occurs if d„ - 6^ v^27cj where / is an integer. Clearly the failure to calculate phase 
a„,„ in such a case is due to the phase shifts 4« and 6m being essentially equal and the 
inter-frame difference function g„„ being identically zero. 

The phase difference between dependent pairs of contingent analytic images, 
such as g„,„ and , is calculated as follows: 



This allows a difference in phase between any two fringe patterns and fk to be 
found, namely {S„-Si(}, with a sign ambiguity as a result of the polarity factor h. It is noted 
that all phase calculations, including phase differences, are calculated modulo 2fr , 
Various methods may be used to extract both the polarity factor h and the phase 
differences In a simplest implementation, from Equation (53), the absolute values 

of the phase difference between two fringe patterns /„ and ft are averaged over the full 
fringe pattern area 5: 



A preferred implementation uses a weighted integral with weighting 'iv(x,yj over 
the fringe pattern area S. The weighting 'i^(x,y) used should be a measure of the estimated 
accuracy of the local estimate of the phase anm- This is related to the AMFM 
characteristics of the fringe pattems f„ix,y) in the sequence of phase-related fringe 
patterns. Thus, the absolute values of the phase difference {S„-Sk) between two fringe 
pattems /^i and/^, weighted averaged over the ftill fringe pattern area S is: 



This allows for a very low (or zero) weighting y^(x,y) to be used where the phase 
a„,n varies quickly, as is the case for discontinuities, or areas where the modulation 
bix,y) IS low- Equation (53) can be rewritten as: 




(54) 




(55) 




(56) 
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By arbitrarily taking the phase difference (<^-4:) for these particular values of n 
and k as positive, from Equation (56) follows an estimated polarity factor he as: 

This estimated polarity factor function he can be used in all further computations, 
because a global sign convention has effectively been set for the phase difference {5„~Sk) 
in the fringe pattern sequence. 

Next, the phase difference S„ - S^^ may be calculated for all possible pairs. The 
minimum reqiiirement is to calculate the phase difference d„ - 6^^ for all sequential pairs, 
which amounts to N values for N fringe patterns f„. By also calculating the phase 
difference 8„ - for non-sequential pairs, additional statistical robustness in the 
estimates is achieved. The number of combinations of two-phase differences 5„~Si^ 
from N phases are: 

N\ 



2{N-2)\ 



(58) 



In the particularly interesting case where N=3, the number of combination is also 
3, which equals the number of sequential differences too. Figs. ISA and 18B show the 
phase shift parameters for the case where N=5 on a polar plane. The points 1-5 are 
thus on a unity circle at points corresponding to exp(ic5rt). The connecting line between 
points n and m corresponds to [exp(j5«)-exp(i5„)]. The length of each connecting line is 
proportional to sin[(<5,„ -<5„)/2], which also appeared in the inter-frame difference 
function g-„„, of Equation (47). The length of each connecting line is a direct indication 
of the reliability of the estimate of the phase difference S„-Sf,, as the length is the factor 
appearing in the magnitude of the inter-frame difference function g^,^ . Fig. 1 8 A shows 
the connecting lines for all sequential pairs, whereas Fig. 18B shows all possible 
connecting lines, and thus inter-frame differences S„—S^., which adds another 5 phase 
differences S^-$^. With the phase differences y„,„ defined as =^ d„-5„ , an 
estimation f of each of the phase differences may be calculated using a weighted 
mean, and compensated by the estimated polarity hg as follows: 
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(59) 



In situations where the firinge patterns /„ contain noise and/or other distortions, 
the best phase differences estimates may come from a two dimensional fit of all the 
points corresponding to exp(/<5„), as shown in Figs. 18A and 18B, to a unit radius circle. 

Without loss of generality, one of the phase shift parameters 6„ may be set to 
zero. For convenience the first phase shift parameter 3^=0 is set. This is consistent 
with the idea that the absolute value of the phase :!(ix,y) is undefined and only phase 
differences y^^ =5^-S^ are really meaningful, The remaining (relative) phase shift 
parameters are calculated as follows: 



These values of the phase shift parameters 5„ are substituted back into the 
general PSA defined by Equation (39). This allows for all the fringe pattern parameters, 
namely the offset a[x,y), the amplitude modulation b{x,y) , and a phase x\X'>y) to be 
calculated, It is noted that the calculated phase x'{>^,y) then includes a common phase 
shift, namely the phase shift S\ . 

However, there may be some residual calibration error due to the original vortex 
operator failing in certain areas. These errors are rectified by re-evaluating the weight 
function w[x,y) based on the present estimates of the offset a{x,y), the amplitude 
modulation b{x,y) , and phase zXx,y), as well as their partial derivatives. This allows 
for the assessment of how well the original fringe patterns, in the form of the inter- frame 
difference functions g^^^ , fulfilled the smoothness requirements of the vortex operator 
V{}. Areas where the phase a„„ varies too quickly, or areas where the modulation 
b(x,y) is too low, the weighting function 'w{x,y) is reduced or set to zero. Once a new 



^=0 



<^2 = + 721 



(60) 
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weight function w{x,y) is formed, the phase calibration process, from Equation (55) 
onwards, is rqpeated, resulting in better estimates for the phase shift parameters S„ . 

The accxxracy of the estimates for tbe phase shift parameters S„ may be 
expressed as the variance of phase difference, as follows: 

wherein 

Care is needed in evaluating Equations (61) and (52), as the phase ;t"(j:,>) is 
periodic and false discontinuities can disrupt the calculations if near -;r or tt 

Other methods for estimating these parameters a„t that specifically incorporate 
the cyclic nature may be used instead of Equation (62), for example: 

= arg ||J^ exp [j |a„„ - a„, |]vt/^ {x, y) dxdy^ (63) 

Fig. 19 shows a preferred embodiment of a automatic phase-shift calibration 
method 600. A first step 610 receives a sequence of fringe patterns /„ . Each of the 
fringe patterns /„ is real valued. An example of one of the fringe patterns f„ from the 
sequence is shown in Fig. 20. Step 620 removes the fringe pattern offset a{x,y) from 
each of the fringe patterns . This is done by forming inter-frame differences g„,„ as 
defined in Equation (40). The inter-frame differences g„„, are also real valued. 

Step 630 forms analytic images g„„, corresponding to the inter- frame differences 
g„„, as defined in Equation (51). This is preferably achieved by placing the inter-frame 
difference function g„„ in the imaginary part of a complex image, and placing the 
imaginary output of the vortex operator V{^„„, > in the real part. Fig. 21 A and 21B show 
the magnitude and phase components respectively of an analytic image g„„ of the 
example fringe pattem/„ shown in Fig. 20. To enable the values of the phase components 
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of the analytic image g„„ to be displayed, the values have been adapted to a range 

127 

[0;255] with 0 representing -/r and 255 representing rj^^ ■ 

Step 140 extracts the phase difference cf„„, -a„^ (modulo 2ti) from two 
dependent analytic images g„„ and . Fig. 22 A shows a phase difference - a^^ 
5 between dependent analytic images g„^„ and of the example fringe pattern f„. Again, 
the values of the phase difference oc^ ~oi^ have been adapted to a range [0;255] with 0 
127 

representing -tt and 255 representing > allowing the values to be displayed. The 

dark areas 710 and 711 have values close to zero on the representation, which 
corresponds to values for close to —n: in the phase difference a^^ -i^,nk> whereas the 

10 light areas 720 and 721 have values close to 255 in the representation, which corresponds 
to values close to ;r in the phase difference (X„^ -c^m* ■ Accordingly, a value of 128 in 
the representation corresponds to a value of 0 in the phase difference a^^^ - a„,i^ . The 
spread of values in the representation is shown in the histogram of Fig. 22B. The peaks 
730 and 731 in the histogram corresponds to phase difference -cCnk values- 

15 Step 150 determines a modulus of the phase difference \oc^„-a^\. Fig. 23 A 

shows the modulus | — \ of the phase difference a„„ - represented in Fig. 22A. 
Again the values have been adapted to a range [0;255] with 0 representing -n and 255 
127 

representing k . As can be seen from Fig. 23A, all the values have a small variation. 

This is confirmed in the histogram shown in Fig. 23 B, wherein the spread of values in the 
20 representation is shown. A single peak 910 in the histogram corresponds to modulus of 
phase difference | t2r„„, ~ a,„^ | values. It is noted Uiat the values have only a small spread 
around the peak value. 

Step 660 calculates the estimated polarity factor he(x^)^ using Equation (37), 
The result is shown in Fig. 24 wherein all regions have values -1 or 1, with the dark areas 
25 corresponding to values of -1 and the light areas corresponding to values of 1 . The phase 
differences f„„ axe calculated using Equation (59) and the relative phase-shift parameters 
= + y„o,-i) estimated in a consistent manner in step 670. 
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The relative phase-shift parameters S„ are then substituted into the general PSA 
defined by Equation (39) and the fringe pattern offset a{x,y), the modulation b{x,y), 
and phase z'i^^y) ^® estimated in step 680. The resultant modulation b{x,y), phase 
X\x,y) and offset a{x,y) of the fringe pattern fnix,y) of the example are shown in Figs. 
5 25A, 25B and 26 respectively. It is noted that Figs. 25A and 26 appear almost uniformly 
white, which corresponds to almost no variation in the values of the modulation b{x,y) 
and offset a [x, y) with respect to coordinates (x^y) for the example shown. 

Step 690 estimates the likely errors In the earlier contingent analytic images g„„ 
§n.k hy means of a direct comparison of the analytic image and the PSA image. 
10 Step 692 determines whether these errors are hi^er than a preset threshold. If the errors 
are smaller than the threshold, the method 600 ends in step 696. If these errors are higher 
than the threshold, the errors are incorporated into the weighting ^{jc,^) in step 694 and 
the method 600 returns to step 670 in an iterative manner, where the new weighting 
'^{x,y) is applied to the estimation of the mean phase-difference f^, , until the errors are 
15 sufjficiently low. In a simplest embodiment, the weighting w{x,y) is set to 1 for co- 
ordinates (x,y) where the comparison is good and 0 where the comparison is poor. 

Fig. 27 shows an example where the comparison of step 690 yields a weighting 
w[x, y) which is 1 .0 where the comparison is good and 0.0 where the comparison is poor. 
It is generally noted that in Equations with integrals over all values of x and y, 
20 such as Equations (41, 54, 55, 59 and 61), the integrals may be replaced by discrete 
summations for discretely sampled images. 

The implementation is applicable to interferometric testing of optical 
components and optical systems. The embodiment may also be used in non-destructive 
testing to study the deformation of objects under stress, and the study of vibrating objects 
25 to identify vibration modes. Often, optical techniques are also used to visuaUse fluid 
flows, and variations in density are displayed in the form of hinge patterns. 

The methods of demodulating a two-dimensional pattern, of estimating an 
orientation angle of a pattern in an image, and of estimating relative phase shifts between 
30 fringe pattern images may be practiced using a conventional general-purpose computer 
system 100, such as that shown in Fig, 28 wherein the processes of Figs. 3, 4, 16 and 19 
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may be implemented as software, such as an application program executing within the 
compulex system 800. In particular, the steps of such methods are effected by instructions 
in the software that are carried out by tiie computer. The software may be stored in a 
computer readable medium, including the storage devices described below^- The software 
is loaded into the computer from the computer readable medium, and then executed by 
the computer. A computer readable medium having such software or computer program 
recorded on it is a computer program product. The use of the computer program product 
in the computer preferably effects an advantageous apparatus for the methods described 
above. 

The computer system 800 comprises a computer module 802, input devices such 
as a keyboard 810 and mouse 812, output devices including a printer 808 and a display 
device 804. 

The computer module 802 typically includes at least one processor unit 814, a 
memory unit 818, for example formed from semiconductor random access memory 
(RAM) and read only memory (ROM), input/output (I/O) interfaces including a video 
interface 822, and an I/O interface 816 for the keyboard 810 and mouse 812. A storage 
device 824 is provided and typically includes a hard disk drive 826 and a floppy disk 
drive 828. A magnetic tape drive (not illustrated) may also be used. A CD-ROM 
drive 820 is typically provided as a non-volatile source of data. The components 814 
to 828 of the computer module 802, typically communicate via an interconnected bus 830 
and in a manner which results in a conventional mode of operation of the computer 
system 100 known to those in the relevant art. Examples of computers on which the 
embodiments can be practised include IBM-PC's and compatibles, Sun Sparcstations or 
alike computer systems evolved therefrom. 

Typically, the appUcation program of the preferred embodiment is resident on 
tlie hard disk drive 826 and read and controlled in its execution by the processor 814. 
Intermediate storage of the program may be accomphshed using the semiconductor 
memory 818, possibly in concert with the hard disk drive 826. In some mstances, the 
application program may be supplied to the user encoded on a CD-ROM or floppy disk 
and read via the corresponding drive 820 or 828, or alternatively may be read by the user 
from a network via a modem device (not illustrated). Still further, the software can also 
be loaded into the computer system 800 from other computer readable medium including 
magnetic tape, a ROM or integrated circuit, a magneto-optical disk, a radio or infra-red 
transmission channel between the computer module 802 and another device, a computer 
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readable card such as a PCMCIA card, and the Internet and Intranets including email 
transmissions and mformation recorded on websites and the like. The foregoing is merely 
exemplary of relevant computer readable mediums. Other computer readable mediums 
may be practiced without departing &om the scope and spirit of the invention. 

The methods may alternatively be implemented in dedicated hardware such as 
one or more integrated circuits performing the functions or sub fimctions. Such dedicated 
hardware may include graphic processors, digital signal processors, or one or more 
microprocessors and associated memories. 

The foregoing describes only some embodiments of the present invention, and 
modifications and/or changes can be made thereto without departing from the scope and 
spirit of the invention, the embodiment being illustrative and not restrictive. 
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